rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-table_A6

##########
#Packages#
##########

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("fastDummies")
library("dplyr")
library("berryFunctions")
library("mediation")
library("sp")
library("pgirmess")
library("fwildclusterboot")
library("geojsonsf")

#Functions

dlshape=function(shploc, shpfile) {
  temp=tempfile()
  download.file(shploc, temp)
  unzip(temp)
  fp <- file.path(temp)
  obj_shp<-rgdal::readOGR(".",shpfile)
  return(obj_shp)}

mean_sd_fun<-function(lm_model) {
  list_empty <-c()
  number_mean<-mean(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  number_sd<-sd(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  list_empty[1]<-number_mean
  list_empty[2]<-number_sd
  names(list_empty)<-c("mean", "sd")
  return(list_empty)
}


week_obs_count<-function(lm_model) {
  list_empty <-c()
  week_no<-length(grep(x = colnames(data.frame(lm_model$model)), pattern = "^week"))
  #Adding 2 because there should be 53 weeks and because I am taking 2 weeks out as in STATA 2 weeks drop
  week_no<-week_no+2
  total_obs<-nobs(lm_model)
  counties<-total_obs/week_no
  list_empty[1]<-week_no
  list_empty[2]<-total_obs
  list_empty[3]<-counties
  names(list_empty)<-c("week_no", "total_obs", "counties")
  return(list_empty)
}


#############################
#Function for calculating SE#
#############################

#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)
data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)


data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)

#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)

#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000
mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)

#Right Wing
data_cov_mob$pct_afd<-as.numeric(data_cov_mob$pct_afd)
mean_pct_afd<-summary(data_cov_mob$pct_afd)[3]
data_cov_mob$pct_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_afd, 1, 0)

#Right Wing Change
data_cov_mob$pct_change_afd<-as.numeric(data_cov_mob$pct_change_afd)
mean_pct_change_afd<-summary(data_cov_mob$pct_change_afd)[4]
data_cov_mob$pct_change_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_change_afd, 1, 0)


data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop", 
                                                        "bonding_tot_pop",
                                                        "pct_afd_dummy",
                                                        "pct_change_afd_dummy"))
data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)

data$post_treat_vereine_tot_pop_dummy<-data$post_treat*data$vereine_tot_pop_dummy
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy
data$post_treat_bridging_tot_pop<-data$post_treat*data$bridging_tot_pop

data$post_treat_bonding_tot_pop_dummy<-data$post_treat*data$bonding_tot_pop_dummy
data$post_treat_pct_afd_dummy<-data$post_treat*data$pct_afd_dummy
data$post_treat_pct_change_afd_dummy<-data$post_treat*data$pct_change_afd_dummy
data$post_treat_pct_change_afd<-data$post_treat*data$pct_change_afd


#Save to separate file
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(nuts2$AGS<-nuts2$AGS)


#For Server too
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')
summary(datas$gender_ratio)


new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))


list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
list_week_dummies_now1<-list_week_dummies[list_week_dummies!= "week_01"]

#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_10"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_53"]


list_land_dummies<-names(new3)
#Re-enable this later if necessary
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]
temp_lag<-10


#Creating the interaction terms for bonding
for (i in list_land_dummies){
  for (j in list_week_dummies_now1)
    datas[[paste("inter_", i, "_", j, sep="")]]<-datas[[i]]*datas[[j]]
}

new<-datas %>% dplyr::select(starts_with("inter_"))
list_land_dummies_star<-names(new)



#BRIDGING
basic_brid<-c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc")
extended_brid <- c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students", "uni_pop")

basic_and_dummies_brid<-list.append(basic_brid, list_dummies, list_week_dummies2)
basic_and_dummies_int_brid<-list.append(basic_brid, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2, list_land_dummies_star)

spec_basic_brid<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_brid, collapse= "+")))
spec_basic_int_brid<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_brid, collapse= "+")))

spec_extended_brid<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_brid, collapse= "+")))
spec_extended_int_brid<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_brid, collapse= "+")))


#Control for Urbanism
#We control for urbanism in a variety of ways
#1) Sample frei stadt and kreis
#2) Above mean population density vs. below mean population density
#In the large models, we already control for pct. student population and no. of universities per capita

datas_frei<-subset(datas, BEZ=="Kreisfreie Stadt" | BEZ=="Stadtkreis")
datas_freix<-subset(datas_frei, week==11)
datas_kreis<-subset(datas, BEZ=="Kreis" | BEZ=="Landkreis")
datas_kreisx<-subset(datas_kreis, week==11)
datas_x<-subset(datas, week==11)

###########################
#FREI STADT LEVEL-ANALYSIS#
###########################

mo1_brid_frei<-lm(spec_basic_int_brid, data=datas_frei)
summary(mo1_brid_frei)
mo1_brid_frei_df<-broom::tidy(mo1_brid_frei)
mo1_brid_frei_df$rse<-coeftest(mo1_brid_frei, vcov = cluster.vcov(mo1_brid_frei, datas_frei$AGS, force_posdef=T))[, 2]
mo1_brid_frei_df$pv_rse<-coeftest(mo1_brid_frei, vcov = cluster.vcov(mo1_brid_frei, datas_frei$AGS, force_posdef=T))[, 4]
mo1_brid_frei_df2<-mo1_brid_frei_df %>% filter(term %in% basic_brid)

#Conley SE
x3<-labels(terms(spec_basic_int_brid))
datas_frei$time<-as.numeric(datas_frei$week)

datas3<-subset(datas_frei, select = list.append(x3, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)


#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_brid, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_brid<-broom::tidy(res)
res_basic_int_brid_df2<-res_basic_int_brid %>% filter(term %in% basic_brid)
res_basic_int_brid_df2<-subset(res_basic_int_brid_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model3<-left_join(mo1_brid_frei_df2, res_basic_int_brid_df2, by = c("term" = "term"))


#####################
#KREIS-LEVEL ANALSYS#
#####################
mo1_brid_kreis<-lm(spec_basic_int_brid, data=datas_kreis)
summary(mo1_brid_kreis)
mo1_brid_kreis_df<-broom::tidy(mo1_brid_kreis)
mo1_brid_kreis_df$rse<-coeftest(mo1_brid_kreis, vcov = cluster.vcov(mo1_brid_kreis, datas_kreis$AGS, force_posdef=T))[, 2]
mo1_brid_kreis_df$pv_rse<-coeftest(mo1_brid_kreis, vcov = cluster.vcov(mo1_brid_kreis, datas_kreis$AGS, force_posdef=T))[, 4]
mo1_brid_kreis_df2<-mo1_brid_kreis_df %>% filter(term %in% basic_brid)


#Conley SE
x3<-labels(terms(spec_basic_int_brid))
datas_kreis$time<-as.numeric(datas_kreis$week)

datas3<-subset(datas_kreis, select = list.append(x3, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas3<-datas3[complete.cases(datas3), ]
datas3_oneweek<-subset(datas3, week_11==1)
datas3_oneweek$AGS<-as.numeric(datas3_oneweek$AGS)

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas3_oneweek<-st_drop_geometry(datas3_oneweek)
nuts_merged<-left_join(nuts2, datas3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)


#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_brid, data=datas3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_brid<-broom::tidy(res)
res_basic_int_brid_df2<-res_basic_int_brid %>% filter(term %in% basic_brid)
res_basic_int_brid_df2<-subset(res_basic_int_brid_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model4<-left_join(mo1_brid_kreis_df2, res_basic_int_brid_df2, by = c("term" = "term"))


#######################
#University per capita#
#######################

mean_uni_pop<-mean(datas$uni_pop, na.rm=TRUE)
datas$above_mean_uni_pop<-ifelse(datas$uni_pop>=mean_uni_pop, 1, 0)
datas_above_mean_uni_pop<-subset(datas, above_mean_uni_pop==1)
datas_below_mean_uni_pop<-subset(datas, above_mean_uni_pop==0)

mo1_brid_above_mean_uni_pop<-lm(spec_basic_int_brid, data=datas_above_mean_uni_pop)
mo1_brid_above_mean_uni_pop_df<-broom::tidy(mo1_brid_above_mean_uni_pop)
mo1_brid_above_mean_uni_pop_df$rse<-coeftest(mo1_brid_above_mean_uni_pop, vcov = cluster.vcov(mo1_brid_above_mean_uni_pop, datas_above_mean_uni_pop$AGS, force_posdef=T))[, 2]
mo1_brid_above_mean_uni_pop_df$pv_rse<-coeftest(mo1_brid_above_mean_uni_pop, vcov = cluster.vcov(mo1_brid_above_mean_uni_pop, datas_above_mean_uni_pop$AGS, force_posdef=T))[, 4]
mo1_brid_above_mean_uni_pop_df2<-mo1_brid_above_mean_uni_pop_df %>% filter(term %in% basic_brid)



mo1_brid_below_mean_uni_pop<-lm(spec_basic_int_brid, data=datas_below_mean_uni_pop)
mo1_brid_below_mean_uni_pop_df<-broom::tidy(mo1_brid_below_mean_uni_pop)
mo1_brid_below_mean_uni_pop_df$rse<-coeftest(mo1_brid_below_mean_uni_pop, vcov = cluster.vcov(mo1_brid_below_mean_uni_pop, datas_below_mean_uni_pop$AGS, force_posdef=T))[, 2]
mo1_brid_below_mean_uni_pop_df$pv_rse<-coeftest(mo1_brid_below_mean_uni_pop, vcov = cluster.vcov(mo1_brid_below_mean_uni_pop, datas_below_mean_uni_pop$AGS, force_posdef=T))[, 4]
mo1_brid_below_mean_uni_pop_df2<-mo1_brid_below_mean_uni_pop_df %>% filter(term %in% basic_brid)


week_no<-length(grep(x = colnames(data.frame(mo1_brid_above_mean_uni_pop$model)), pattern = "^week")) + 2
total_obs<-nobs(mo1_brid_above_mean_uni_pop)
counties<-total_obs/week_no
week_no<-length(grep(x = colnames(data.frame(mo1_brid_below_mean_uni_pop$model)), pattern = "^week")) + 2
total_obs<-nobs(mo1_brid_below_mean_uni_pop)
counties<-total_obs/week_no


#######################
#University Below Mean#
#######################

mo1_brid_below_mean_uni_pop<-lm(spec_basic_int_brid, data=datas_below_mean_uni_pop)
summary(mo1_brid_below_mean_uni_pop)
mo1_brid_below_mean_df<-broom::tidy(mo1_brid_below_mean_uni_pop)
mo1_brid_below_mean_df$rse<-coeftest(mo1_brid_below_mean_uni_pop, vcov = cluster.vcov(mo1_brid_below_mean_uni_pop, datas_below_mean_uni_pop$AGS, force_posdef=T))[, 2]
mo1_brid_below_mean_df$pv_rse<-coeftest(mo1_brid_below_mean_uni_pop, vcov = cluster.vcov(mo1_brid_below_mean_uni_pop, datas_below_mean_uni_pop$AGS, force_posdef=T))[, 4]
mo1_brid_below_mean_df2<-mo1_brid_below_mean_df %>% filter(term %in% basic_brid)


x3<-labels(terms(spec_basic_int_brid))
datas_below_mean_uni_pop$time<-as.numeric(datas_below_mean_uni_pop$week)

datas_below_mean_uni_pop3<-subset(datas_below_mean_uni_pop, select = list.append(x3, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas_below_mean_uni_pop3<-datas_below_mean_uni_pop3[complete.cases(datas_below_mean_uni_pop3), ]
datas_below_mean_uni_pop3_oneweek<-subset(datas_below_mean_uni_pop3, week_11==1)
datas_below_mean_uni_pop3_oneweek$AGS<-as.numeric(datas_below_mean_uni_pop3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas_below_mean_uni_pop3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas_below_mean_uni_pop3_oneweek<-st_drop_geometry(datas_below_mean_uni_pop3_oneweek)
nuts_merged<-left_join(nuts2, datas_below_mean_uni_pop3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)



#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas_below_mean_uni_pop3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_brid, data=datas_below_mean_uni_pop3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_brid<-broom::tidy(res)
res_basic_int_brid_df2<-res_basic_int_brid %>% filter(term %in% basic_brid)
res_basic_int_brid_df2<-subset(res_basic_int_brid_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')

model16<-left_join(mo1_brid_below_mean_df2, res_basic_int_brid_df2, by = c("term" = "term"))


#######################
#University Below Mean#
#######################


mo1_brid_above_mean_uni_pop<-lm(spec_basic_int_brid, data=datas_above_mean_uni_pop)
summary(mo1_brid_above_mean_uni_pop)
mo1_brid_above_mean_uni_pop_df<-broom::tidy(mo1_brid_above_mean_uni_pop)
mo1_brid_above_mean_uni_pop_df$rse<-coeftest(mo1_brid_above_mean_uni_pop, vcov = cluster.vcov(mo1_brid_above_mean_uni_pop, datas_above_mean_uni_pop$AGS, force_posdef=T))[, 2]
mo1_brid_above_mean_uni_pop_df$pv_rse<-coeftest(mo1_brid_above_mean_uni_pop, vcov = cluster.vcov(mo1_brid_above_mean_uni_pop, datas_above_mean_uni_pop$AGS, force_posdef=T))[, 4]
mo1_brid_above_mean_uni_pop_df2<-mo1_brid_above_mean_uni_pop_df %>% filter(term %in% basic_brid)


#######################
#University above Mean#
#######################

mo1_brid_above_mean_uni_pop<-lm(spec_basic_int_brid, data=datas_above_mean_uni_pop)
summary(mo1_brid_above_mean_uni_pop)
mo1_brid_above_mean_df<-broom::tidy(mo1_brid_above_mean_uni_pop)
mo1_brid_above_mean_df$rse<-coeftest(mo1_brid_above_mean_uni_pop, vcov = cluster.vcov(mo1_brid_above_mean_uni_pop, datas_above_mean_uni_pop$AGS, force_posdef=T))[, 2]
mo1_brid_above_mean_df$pv_rse<-coeftest(mo1_brid_above_mean_uni_pop, vcov = cluster.vcov(mo1_brid_above_mean_uni_pop, datas_above_mean_uni_pop$AGS, force_posdef=T))[, 4]
mo1_brid_above_mean_df2<-mo1_brid_above_mean_df %>% filter(term %in% basic_brid)


x3<-labels(terms(spec_basic_int_brid))
datas_above_mean_uni_pop$time<-as.numeric(datas_above_mean_uni_pop$week)

datas_above_mean_uni_pop3<-subset(datas_above_mean_uni_pop, select = list.append(x3, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y", "time", "AGS"))
datas_above_mean_uni_pop3<-datas_above_mean_uni_pop3[complete.cases(datas_above_mean_uni_pop3), ]
datas_above_mean_uni_pop3_oneweek<-subset(datas_above_mean_uni_pop3, week_11==1)
datas_above_mean_uni_pop3_oneweek$AGS<-as.numeric(datas_above_mean_uni_pop3_oneweek$AGS)

#Conley SE

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(datas_above_mean_uni_pop3$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)


#Merging with spatial dataframe
datas_above_mean_uni_pop3_oneweek<-st_drop_geometry(datas_above_mean_uni_pop3_oneweek)
nuts_merged<-left_join(nuts2, datas_above_mean_uni_pop3_oneweek, by = c("AGS"="AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)



#Calculating the distance the minimises Moran's I stat
dist_cut_brid<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

dm_bri <- dist_mat(datas_above_mean_uni_pop3, dist_cutoff = dist_cut_brid, 
                   lat = "POINT_Y", lon = "POINT_X", unit = "AGS", time = "time",
                   dist_round = TRUE)
res<-conleyreg(spec_basic_int_brid, data=datas_above_mean_uni_pop3, dist_cutoff=dist_cut_brid, unit = "AGS", time = "time", 
               dist_mat=dm_bri, lag_cutoff=temp_lag)

res_basic_int_brid<-broom::tidy(res)
res_basic_int_brid_df2<-res_basic_int_brid %>% filter(term %in% basic_brid)
res_basic_int_brid_df2<-subset(res_basic_int_brid_df2, select = c("term", "std.error", 'p.value'))
names(res_basic_int_brid_df2)<-c("term", "con_std_error", 'con_p_value')


model17<-left_join(mo1_brid_above_mean_df2, res_basic_int_brid_df2, by = c("term" = "term"))


#Arranging models

model3$order<-NA
model4$order<-NA
model16$order<-NA
model17$order<-NA


model3$order[model3$term=="bridging_tot_pop_dummy"]<-1
model4$order[model4$term=="bridging_tot_pop_dummy"]<-1
model16$order[model16$term=="bridging_tot_pop_dummy"]<-1
model17$order[model17$term=="bridging_tot_pop_dummy"]<-1

model3$order[model3$term=="post_treat"]<-2
model4$order[model4$term=="post_treat"]<-2
model16$order[model16$term=="post_treat"]<-2
model17$order[model17$term=="post_treat"]<-2

model3$order[model3$term=="post_treat_bridging_tot_pop_dummy"]<-3
model4$order[model4$term=="post_treat_bridging_tot_pop_dummy"]<-3
model16$order[model16$term=="post_treat_bridging_tot_pop_dummy"]<-3
model17$order[model17$term=="post_treat_bridging_tot_pop_dummy"]<-3


model3$estimate<-round(as.numeric(model3$estimate), 3)
model4$estimate<-round(as.numeric(model4$estimate), 3)
model16$estimate<-round(as.numeric(model16$estimate), 3)
model17$estimate<-round(as.numeric(model17$estimate), 3)

model3$se<-paste("(", round(as.numeric(model3$std.error),3), ")",
                 ifelse(as.numeric(model3$p.value)<0.001,"***",
                        ifelse(as.numeric(model3$p.value<0.01),"**",
                               ifelse(as.numeric(model3$p.value)<0.05,"*",
                                      ifelse(as.numeric(model3$p.value)<0.1,"+",
                                      "")))), sep = "")

model4$se<-paste("(", round(as.numeric(model4$std.error),3), ")",
                 ifelse(as.numeric(model4$p.value)<0.001,"***",
                        ifelse(as.numeric(model4$p.value<0.01),"**",
                               ifelse(as.numeric(model4$p.value)<0.05,"*",
                                      ifelse(as.numeric(model4$p.value)<0.1,"+",
                                             "")))), sep = "")


model16$se<-paste("(", round(as.numeric(model16$std.error),3), ")",
                  ifelse(as.numeric(model16$p.value)<0.001,"***",
                         ifelse(as.numeric(model16$p.value<0.01),"**",
                                ifelse(as.numeric(model16$p.value)<0.05,"*",
                                       ifelse(as.numeric(model16$p.value)<0.1,"+",
                                              "")))), sep = "")

model17$se<-paste("(", round(as.numeric(model17$std.error),3), ")",
                  ifelse(as.numeric(model17$p.value)<0.001,"***",
                         ifelse(as.numeric(model17$p.value<0.01),"**",
                                ifelse(as.numeric(model17$p.value)<0.05,"*",
                                       ifelse(as.numeric(model17$p.value)<0.1,"+",
                                              "")))), sep = "")


model3$std_er_star<-paste("(", round(as.numeric(model3$rse),3), ")",
                          ifelse(as.numeric(model3$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model3$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model3$p.value)<0.05,"*",
                                               ifelse(as.numeric(model3$p.value)<0.1,"+",
                                                      "")))), sep = "")

model4$std_er_star<-paste("(", round(as.numeric(model4$rse),3), ")",
                          ifelse(as.numeric(model4$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model4$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model4$p.value)<0.05,"*",
                                               ifelse(as.numeric(model4$p.value)<0.1,"+",
                                                      "")))), sep = "")


model16$std_er_star<-paste("(", round(as.numeric(model16$rse),3), ")",
                           ifelse(as.numeric(model16$pv_rse)<0.001,"***",
                                  ifelse(as.numeric(model16$pv_rse<0.01),"**",
                                         ifelse(as.numeric(model16$p.value)<0.05,"*",
                                                ifelse(as.numeric(model16$p.value)<0.1,"+",
                                                       "")))), sep = "")

model17$std_er_star<-paste("(", round(as.numeric(model17$rse),3), ")",
                           ifelse(as.numeric(model17$pv_rse)<0.001,"***",
                                  ifelse(as.numeric(model17$pv_rse<0.01),"**",
                                         ifelse(as.numeric(model17$p.value)<0.05,"*",
                                                ifelse(as.numeric(model17$p.value)<0.1,"+",
                                                       "")))), sep = "")


model3$con_std_er_star<-paste("[", round(as.numeric(model3$con_std_error),3), "]",
                              ifelse(as.numeric(model3$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model3$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model3$p.value)<0.05,"*",
                                                   ifelse(as.numeric(model3$p.value)<0.1,"+",
                                                          "")))), sep = "")


model4$con_std_er_star<-paste("[", round(as.numeric(model4$con_std_error),3), "]",
                              ifelse(as.numeric(model4$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model4$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model4$p.value)<0.05,"*",
                                                   ifelse(as.numeric(model4$p.value)<0.1,"+",
                                                          "")))), sep = "")


model16$con_std_er_star<-paste("[", round(as.numeric(model16$con_std_error),3), "]",
                              ifelse(as.numeric(model16$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model16$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model16$p.value)<0.05,"*",
                                                   ifelse(as.numeric(model16$p.value)<0.1,"+",
                                                          "")))), sep = "")

model17$con_std_er_star<-paste("[", round(as.numeric(model17$con_std_error),3), "]",
                              ifelse(as.numeric(model17$con_p_value)<0.001,"***",
                                     ifelse(as.numeric(model17$con_p_value<0.01),"**",
                                            ifelse(as.numeric(model17$p.value)<0.05,"*",
                                                   ifelse(as.numeric(model17$p.value)<0.1,"+",
                                                          "")))), sep = "")



model3_b<-subset(model3, select = c("order", "term", "estimate", "se", "std_er_star", "con_std_er_star"))
samples<-c(20, "sample", "Free Cities (Kreisfreie Stadt)", "", "", "", "")
means<-c(21, "mean_mob", round(mean_sd_fun(mo1_brid_frei)[1], 3), "", "", "", "")
sds<-c(22, "sd_mob", round(mean_sd_fun(mo1_brid_frei)[2], 3), "", "", "", "")
no_cross_sections<-c(23, "cross_sections", round(week_obs_count(mo1_brid_frei)[3], 3), "", "", "", "")
no_time_units<-c(24, "time_units", round(week_obs_count(mo1_brid_frei)[1], 3), "", "", "", "")
county_fe<-c(25, "county_fe", "Yes", "", "", "", "")
week_fe<-c(26, "week_fe", "Yes", "", "", "", "")
weekxland_fe<-c(27, "weekxland_fe", "Yes", "", "", "", "")
adj_rsq<-c(28, "adj_rsq", round(summary(mo1_brid_frei)$adj.r.squared, 3), "", "", "", "")
obs<-c(29, "obs", nobs(mo1_brid_frei), "", "", "", "")

model3_b<-rbind(model3_b, samples, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
model3_b$model<-"Model1"
model3_b$order<-as.numeric(model3_b$order)
names(model3_b)<-paste0(names(model3_b), "_1")



model4_b<-subset(model4, select = c("order", "term", "estimate", "se", "std_er_star", "con_std_er_star"))
samples<-c(20, "sample", "Counties (Kreis)", "", "", "", "")
means<-c(21, "mean_mob", round(mean_sd_fun(mo1_brid_kreis)[1], 3), "", "", "", "")
sds<-c(22, "sd_mob", round(mean_sd_fun(mo1_brid_kreis)[2], 3), "", "", "", "")
no_cross_sections<-c(23, "cross_sections", round(week_obs_count(mo1_brid_kreis)[3], 3), "", "", "", "")
no_time_units<-c(24, "time_units", round(week_obs_count(mo1_brid_kreis)[1], 3), "", "", "", "")
county_fe<-c(25, "county_fe", "Yes", "", "", "", "")
week_fe<-c(26, "week_fe", "Yes", "", "", "", "")
weekxland_fe<-c(27, "weekxland_fe", "Yes", "", "", "", "")
adj_rsq<-c(28, "adj_rsq", round(summary(mo1_brid_kreis)$adj.r.squared, 3), "", "", "", "")
obs<-c(29, "obs", nobs(mo1_brid_kreis), "", "", "", "")

model4_b<-rbind(model4_b, samples, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
model4_b$model<-"Model1"
model4_b$order<-as.numeric(model4_b$order)
names(model4_b)<-paste0(names(model4_b), "_2")





model16_b<-subset(model16, select = c("order", "term", "estimate", "se", "std_er_star", "con_std_er_star"))
samples<-c(20, "sample", "Above Mean Uni-Pop. Ratio", "", "", "", "")
means<-c(21, "mean_mob", round(mean_sd_fun(mo1_brid_above_mean_uni_pop)[1], 3), "", "", "", "")
sds<-c(22, "sd_mob", round(mean_sd_fun(mo1_brid_above_mean_uni_pop)[2], 3), "", "", "", "")
no_cross_sections<-c(23, "cross_sections", round(week_obs_count(mo1_brid_above_mean_uni_pop)[3], 3), "", "", "", "")
no_time_units<-c(24, "time_units", round(week_obs_count(mo1_brid_above_mean_uni_pop)[1], 3), "", "", "", "")
county_fe<-c(25, "county_fe", "Yes", "", "", "", "")
week_fe<-c(26, "week_fe", "Yes", "", "", "", "")
weekxland_fe<-c(27, "weekxland_fe", "Yes", "", "", "", "")
adj_rsq<-c(28, "adj_rsq", round(summary(mo1_brid_above_mean_uni_pop)$adj.r.squared, 3), "", "", "", "")
obs<-c(29, "obs", nobs(mo1_brid_above_mean_uni_pop), "", "", "", "")

model16_b<-rbind(model16_b, samples, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
model16_b$model<-"Model1"
model16_b$order<-as.numeric(model16_b$order)
names(model16_b)<-paste0(names(model16_b), "_5")



model17_b<-subset(model17, select = c("order", "term", "estimate", "se", "std_er_star", "con_std_er_star"))
samples<-c(20, "sample", "Below Mean Uni-Pop. Ratio", "", "", "", "")
means<-c(21, "mean_mob", round(mean_sd_fun(mo1_brid_below_mean_uni_pop)[1], 3), "", "", "", "")
sds<-c(22, "sd_mob", round(mean_sd_fun(mo1_brid_below_mean_uni_pop)[2], 3), "", "", "", "")
no_cross_sections<-c(23, "cross_sections", round(week_obs_count(mo1_brid_below_mean_uni_pop)[3], 3), "", "", "", "")
no_time_units<-c(24, "time_units", round(week_obs_count(mo1_brid_below_mean_uni_pop)[1], 3), "", "", "", "")
county_fe<-c(25, "county_fe", "Yes", "", "", "", "")
week_fe<-c(26, "week_fe", "Yes", "", "", "", "")
weekxland_fe<-c(27, "weekxland_fe", "Yes", "", "", "", "")
adj_rsq<-c(28, "adj_rsq", round(summary(mo1_brid_below_mean_uni_pop)$adj.r.squared, 3), "", "", "", "")
obs<-c(29, "obs", nobs(mo1_brid_below_mean_uni_pop), "", "", "", "")

model17_b<-rbind(model17_b, samples, means, sds, no_cross_sections, no_time_units, county_fe,
                 week_fe, weekxland_fe, adj_rsq, obs)
model17_b$model<-"Model1"
model17_b$order<-as.numeric(model17_b$order)
names(model17_b)<-paste0(names(model17_b), "_6")





#Step1: Collecting the models
list_order<-c(model3_b$order_1, model4_b$order_2, 
              model16_b$order_5, model17_b$order_6)
#Step2: Ordering rows
list_order<-unique(sort(list_order))
list_order_df<-data.frame(list_order)
names(list_order_df)<-"order"
list_order_df$order<-as.numeric(list_order_df$order)


new_df2<-left_join(list_order_df, as.data.frame(model3_b), by = c("order"="order_1"))


new_df3<-left_join(new_df2, as.data.frame(model4_b), by = c("order"="order_2"))
new_df3$model<-new_df3$model_1
new_df3$model<-ifelse(is.na(new_df3$model), new_df3$model_2, new_df3$model)
new_df3$term<-new_df3$term_1
new_df3$term<-ifelse(is.na(new_df3$term_1), new_df3$term_2, new_df3$term)
new_df3<-subset(new_df3, select = -c(term_1, model_1, term_2, model_2))



new_df6<-left_join(new_df3, as.data.frame(model16_b), by = c("order"="order_5"))
new_df6$model<-ifelse(is.na(new_df6$model), new_df6$model_5, new_df6$model)
new_df6$term<-ifelse(is.na(new_df6$term), new_df6$term_5, new_df6$term)
new_df6<-subset(new_df6, select = -c(term_5, model_5))

new_df7<-left_join(new_df6, as.data.frame(model17_b), by = c("order"="order_6"))
new_df7$model<-ifelse(is.na(new_df7$model), new_df7$model_6, new_df7$model)
new_df7$term<-ifelse(is.na(new_df7$term), new_df7$term_6, new_df7$term)
new_df7<-subset(new_df7, select = -c(term_6, model_6))



new_df7$formal<-NA
new_df7$formal[new_df7$term=="bridging_tot_pop_dummy"]<-"Bridging"
new_df7$formal[new_df7$term=="post_treat"]<-"Post Week 11 Dummy"
new_df7$formal[new_df7$term=="post_treat_bridging_tot_pop_dummy"]<-"Post Week 11 x Bridging"
new_df7$formal[new_df7$term=="samples"]<-"Sample"
new_df7$formal[new_df7$term=="mean_mob"]<-"Mean Mobility"
new_df7$formal[new_df7$term=="sd_mob"]<-"SD Mobility"
new_df7$formal[new_df7$term=="cross_sections"]<-"No. Cross-Sections"
new_df7$formal[new_df7$term=="time_units"]<-"No. Time Units"
new_df7$formal[new_df7$term=="county_fe"]<-"County FE"
new_df7$formal[new_df7$term=="week_fe"]<-"Week FE"
new_df7$formal[new_df7$term=="weekxland_fe"]<-"Week X Land FE"
new_df7$formal[new_df7$term=="adj_rsq"]<-"Adj. R sq."
new_df7$formal[new_df7$term=="obs"]<-"Observations"

new_df7<-subset(new_df7, order!=20)


#Without pop den

new_df9<-subset(new_df7, select = c(formal, estimate_1, se_1, std_er_star_1, con_std_er_star_1,
                                    estimate_2, se_2, std_er_star_2, con_std_er_star_2,
                                    estimate_5, se_5, std_er_star_5, con_std_er_star_5,
                                    estimate_6, se_6, std_er_star_6, con_std_er_star_6))

#Step3 Order columns
final1<-new_df9
new_row<-c("Variable", rep(c("Estimate", "SE", "Robust SE", "Conley SE"), 4))
final2<-rbind(new_row, final1)

# Get the number of columns in new_df8
num_columns <- ncol(new_df9)+1
# Create an alignment string with 'l's for each column
alignment_string <- paste(rep("l", num_columns), collapse = "")
nchar(alignment_string)

#Step7 centering coefficients
object_latex<-xtable(final2, type = "latex", 
                     caption = "Results",
                     #Note there are X columns to be aligned below.
                     #This is because of the index column.
                     align = alignment_string)


comment <- list(pos = list(0), command = NULL)
comment$pos[[1]] <- c(nrow(object_latex))
comment$command <- c(paste("\\bottomrule\n",
                           "\\multicolumn{", num_columns-1, "}{l} {\\parbox[t]{37cm}\n",
                           "{\\footnotesize \\textit{Note}: 
This table shows the effect of networks on mobility in the post-lockdown period. 
The dependent variable is mobility as measured by the government. 
The interactions between the period after week 11 and bridgingg networks. These take the value of 1 for counties with 
density of these individual networks above the mean. Similarly, the post-lockdown period is marked with 1.
The standard errors are clustered at a county level in round brackets and 
adjusted for spatial (bridging - 538 km) and serial (10 weeks) correlation in the square brackets. 
Significant at +10\\%, ∗5\\%, ∗1\\%, and ∗0.1\\%.}} \\\\ \n", sep = ""))

#Step8 Printing latex tex
tableLines<-print(object_latex,
                  #The following line controls where the lines are drawn
                  hline.after=c(-1, 1, 4),
                  floating = FALSE,
                  file = "./paper/tables/table_A6.tex", 
                  #caption.placement = "top",
                  add.to.row = comment,
                  booktabs = TRUE,
                  include.colnames=F,
                  include.rownames=FALSE,
                  sanitize.text.function=function(x){x})

multicolumns <- "& \\\\multicolumn{4}{c}{M1: Free Cities (Kreisfreie Stadt)} & 
\\\\multicolumn{4}{c}{M2: Counties (Kreis)} & 
\\\\multicolumn{4}{c}{M3: Above Mean Uni-Pop. Ratio} &
\\\\multicolumn{4}{c}{M4: Below Mean Uni-Pop. Ratio}\\\\\\\\"
clines<-"\\\\cmidrule(lr){2-5} \\\\cmidrule(lr){6-9} \\\\cmidrule(lr){10-13} \\\\cmidrule(lr){14-17}
\\\\cmidrule(lr){18-21} \\\\cmidrule(lr){22-25}\\\\\\\\"
tableLines <- sub ("\\\\toprule\\n", paste0 ("\\\\toprule\n", multicolumns, clines, "\n"), tableLines) ## booktabs = TRUE
#tableLines <- sub ("\\\\hline\\n",   paste0 ("\\\\hline\n",   multicolumns, "\n"), tableLines) ## booktabs = FALSE
writeLines(tableLines, con = "./paper/tables/table_A6.tex", )



